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We present a computational method to solve the magnetohydrodynamic equations in spherical 
geometry. The technique is fully nonlinear and wholly spectral, and uses an expansion basis that is 
adapted to the geometry: Chandrasekhar-Kendall vector eigenfunctions of the curl. The resulting 
lower spatial resolution is somewhat offset by being able to build all the boundary conditions into 
each of the orthogonal expansion functions and by the disappearance of any difficulties caused 
by singularities at the center of the sphere. The results reported here are for mechanically and 
magnetically isolated spheres, although different boundary conditions could be studied by adapting 
the same method. The intent is to be able to study the nonlinear dynamical evolution of those 
aspects that are peculiar to the spherical geometry at only moderate Reynolds numbers. The code 
is parallelized, and will preserve to high accuracy the ideal magnetohydrodynamic (MHD) invariants 
of the system (global energy, magnetic helicity, cross helicity) . Examples of results for selective decay 
and mechanically-driven dynamo simulations are discussed. In the dynamo cases, spontaneous flips 
of the dipole orientation are observed. 

PACS numbers: 47.ll.-j; 47.11. Kb; 91.25.Cw; 95.30.Qd 



I. INTRODUCTION 

Magnetohydrodynamic "dynamo" processes are those 
in which the motions of an electrically conducting fluid 
amplify and maintain a finite magnetic field, starting 
from an arbitrarily small one. They have long been of 
interest for geophysics and astrophysics Q, S S S IE 111 j 
and have more recently become of interest with regard to 
laboratory attempts to generate dynamo magnetic fields 
in hquid metals IlHSIlElIllIl- The relevant theoreti- 
cal and computational literature is vast, and extensive re- 
views have recently been given (see e.g. Rcfs. P. IT3 . ll^ ). 

In our own work, we have lately been studying dynamo 
processes numerically, using turbulent three-dimensional 
magnetohydrodynamic (hereafter, "MHD") codes of the 
familiar Orszag-Patterson pseudospectral variety 0, 0, 
Ir^ ITM H^ . Such codes treat homogeneous turbulence 
efficiently, but are mainly useful in situations involving 
spatially periodic boundary conditions. Particularly for 
the case of planetary dynamos and laboratory experi- 
ments, restrictions to periodic boundary conditions arc a 
severe limitation. Essential ingredients, such as rotation, 
global angular momentum, and the interfaces between 
conducting and non-conducting regions are not readily 
included. It is to be expected that all of these play a 
role in the physical situations of interest, and give rise to 
qualitatively new physical processes not accessible with 
periodic boundaries. 

This present paper represents our attempt to begin 
a study of these processes by introducing, and display- 
ing some results from, a computational method that is 
adapted to the geometry of isolated spheres. Our goal 
is not to reach realistic geophysical parameter regimes 
(out of the question for the foreseeable future, in any 
case) , but rather to isolate and study those new physical 
processes that appear in this geometrically more realis- 



tic setting. Our results are not to be compared with the 
ambitious geo-dynamo computation of Glatzmaier and 
Roberts (e.g., Refs. [2,21]) who include far more physical 
effects and degrees of freedom than we have at this point. 
We will consider it not to be a limitation on our work if we 
are unable to reach realistic Reynolds numbers, so long 
as the spectra of the dynamo processes we do identify 
are well-resolved by the number of expansion functions 
we arc able to retain. We think of these efforts as studies 
of dynamo behaviors exhibited by the MHD equations 
in spherical geometry without regard to their present 
applicability to geophysical or laboratory situations. It 
appears not to be necessary to reach presently-existing 
ranges of Reynolds numbers or magnetic Prandtl num- 
bers in order to see interesting spontaneously-generated 
magnetic field behavior. 

The numerical technique is entirely spectral, using an 
expansion basis that is specifically adapted to spherical 
geometry: the Chandrasekhar-Kendall (hereafter, "C- 
K") vector eigenfunctions of the curl [20 l |2 L .22 . .23. |23| 
(by construction, C-K functions are als o eig enfunctions of 
the Helmholtz wave equation; see e.g. [23 )• These func- 
tions form a complete orthogonal set under the bound- 
ary conditions we will choose. They can be normal- 
ized. Their conipleteness has been shown by Yoshida 
for the cylinder [23 and by Cantarella et al |2J| for the 
sphere. The cylindrical version was used some years ago 
for studying processes believed to be involved in the dis- 
ruptions of fusion confinement devices |2^ . l2^l23 . l29ll30l| . 
Their advantages lie in their natural geometrical relation 
to the specific geometries in which they are employed (all 
the boundary conditions can be built into each expansion 
function) and in certain desirable numerical properties, 
which include the following. The MHD equations involve 
several solenoidal fields of which several curls are taken. 
Taking the curl of a C-K function simply gives the func- 
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tion back again, times a multiplicative constant. This 
means that no numerical spatial differentiations are re- 
quired, with their attendant complication of the expres- 
sions and loss of accuracy. The code preserves the ideal 
MHD invariants very accurately over long times (total 
energy, magnetic helicity, cross helicity), which is one 
of the few accuracy checks available to a strongly non- 
linear code for which non-trivial analytical solutions are 
scarce. In addition, the ideal invariants are readily ex- 
hibited as simple algebraic quadratic sums involving the 
expansion coefficients, so no numerical integrations are 
required to evaluate them. Going along with these ad- 
vantages is a severe disadvantage associated with the lack 
of a fast transform for the spherical Bessel functions in- 
volved, which means that the nonlinear terms lead to 
convolution sums which grow rapidly with the number of 
expansion functions retained in the Galerkin approxima- 
tion. This limits us to mechanical and magnetic Reynolds 
numbers of the order of a very few thousand, and makes 
it unlikely that we can ever reach planetary parameters 
without modeling the small scales of the fluid motions. 

But the absence of any coordinate singularities (e.g., 
r = in spherical polar coordinates) to worry about any- 
where is a considerable advantage, as we noted previously 
in solving the two-dimensional Navier-Stokes equation in- 
side a circle with non- ideal boundaries |^ [s^ • The code 
also can be readily parallelized and there are no potential 
aliasing problems. 

The situation we wish to study is that of an electri- 
cally conducting fluid inside a rigid spherical boundary 
that isolates the fluid, mechanically and magnetically, 
from everything outside it. The boundary is regarded as 
a spherical shell that is perfectly conducting (so magnetic 
fields cannot penetrate the region outside) and is coated 
on the inside with a thin layer of insulating dielectric so 
that the electrical current density cannot penetrate the 
shell. The shell is regarded as mechanically impenetra- 
ble, so that the normal component of the fluid veloc- 
ity vanishes there. We do not employ no-slip boundary 
conditions on the velocity field, but rather choose the 
vanishing of the normal component of the vorticity at 
the spherical boundary as the second boundary condi- 
tion on the velocity field. This condition is implied by, 
but does not imply, no-slip boundary conditions on the 
velocity field, and thereby avoids certain mathematical 
procedures that sometimes attend the attempts at impos- 
ing no-slip boundary conditions. These four boundary 
conditions for the fields in the external walls, plus reg- 
ularity of each component of the velocity and magnetic 
fields at the origin, might be thought of as a total of ten 
boundary conditions for the system. The above men- 
tioned boundary condition for the vorticity in the wall 
is easy to implement in the present formulation. Other 
boundary conditions (e.g. no-slip boundary conditions) 
can be studied with our method using a penalty method 
close to the walls (as done for example in Rcf. [s^l) and 
will be considered in the future. However, we want to 
point out that the present election of boundary condi- 



tions also avoids some problems present in hydrodynamic 
incompressible fiows when no-slip conditions are imposed 
(see e.g. [sj for a detailed discussion). Since this 
topic is beyond the aim of our present study, we will use 
in the following these simple boundary conditions and 
leave the study of other options for future work. 

An outline of the paper follows. Section HTl introduces 
the expansion basis and shows how the nonlinear partial 
differential equations can be reduced to a set of ordi- 
nary differential equations, first order in the time, for 
the expansion coefficients. We include mention of other 
possible boundary conditions that it is intended to intro- 
duce in the future, then settle on the one just described 
for these runs. The numerical method is discussed in 
Section UTTI Sections IIVI and Ivl displav our first applica- 
tions of the code to some simple MHD problems that are 
considered interesting and that are peculiar to spherical 
geometry. We simultaneously explore the possibilities for 
some relative new flow visualization diagnostics when ex- 
hibiting our results [s^; these are described when they 
appear. Section VI summarizes the results and discusses 
possible future applications of the method. 

II. THE SPECTRAL DECOMPOSITION 

C-K functions are constructed from a solution of the 
scalar Helmholtz equation: 

(v2 + a2)^ = o, (1) 

where we are referring to spherical polar coordinates 
(r, 6, (j)) and A is an eigenvalue that will eventually be de- 
termined by boundary conditions. Vector cigcnfunctions 
of the curl appropriate to spherical geometry may be con- 
structed according to the following recipe from each so- 
lution to Eq. CQ: 

J = AV X r?/- + V X (V X ri/O . (2) 

From this, it may readily be verifled that V x J = AJ. 
Thus any single J is a "force-free" or "Bernoulli" flcld, 
though the sum of two or more of them is not. The 
relevant scalar i}} is 

V'gim = Cqljl{\\ql\r)Yiyn{e,(l)). (3) 

Here, Cqi is a normalization constant, ji is a spherical 
Bessel function of order I, and Yim is the normalized 
spherical harmonic expressed in terms of the polar angle 
9 and the azimuthal (longitudinal) angle (j). The number 
Cqi is chosen to make the volume integral of 3qim ■ Jg/m 
over the computational domain equal to unity (the aster- 
isk denotes complex conjugate). 

The integer / is 1, 2, 3, . . . and m runs in integer steps 
from —I to /. An infinite sequence of values of A^;, la- 
beled by q = 1, 2, 3, . . . and q = —1, —2, —3, . . . will be 
determined by a radial boundary condition momentarily 
to be invoked on J for a given I. 
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The boundary conditions that will be invoked, consis- 
tently with the physical description of the region and its 
boundaries that were given in Sec. prefer to the magnetic 
field B and the velocity field v, both solenoidal, and their 
curls. The curl of v is a;, the vorticity, and the curl of 
B is j, the electric current density, in the dimensionless 
units that remain to be described. We shall demand that 
the normal (radial) components of all four vector fields 
vanish at the boundary r = R. Since these four fields 
are to be expanded in the J's, making the radial compo- 
nent of each J vanish at r = i? will guarantee that the 
boundary conditions will be satisfied for any superposi- 
tion of them. This vanishing of the normal components 
of the J determines the allowed values of Xqi, positive 
and negative, by the locations of the zeros of the spheri- 
cal Bcsscl functions. Then the explicit expression for the 
normalization constant Cqi is 

Cqi = Wui+iil^RT' [1(1 + m'] '^'^ ■ (4) 

Positive or negative q is to be associated with positive or 
negative Ag/ according to A_g./ ~ For a variety of 

radial boundary conditions, including the one just cho- 
sen, the 3qim functions corresponding to differing indices 
q, I, or m can be shown to be orthogonal. Boundary con- 
ditions at two different radii can be imposed if spherical 
Neumann functions are also permitted in ij:. 

Given the completeness of the C-K functions so de- 
fined, they are then appropriate for expanding various 
vector fields of incompressible MHD: the magnetic field, 
fiuid velocity, vorticity, electric current density, and vec- 
tor potential in the Coulomb gauge. These are consid- 
ered to be the necessary set of quantities to be expanded 
when studying spherical MHD and dynamos in the class 
of problems studied here. The ones of slowest spatial 
variation (smaller low I and |m|) contain the dipole 
and low-order multipole components of the fields and pro- 
vide a natural ordering of the spectral contributions from 
various spatial scales, starting with those of the largest 
scales. 

The MHD equations to be solved are, in familiar di- 
mensionless ("Alfvcnic") units, common in MHD turbu- 
lence theory, an equation of motion for the fiuid velocity 

V, 

— = V X X B- V f 7^ + y j +i^V2v + f, (5) 

and the induction equation for advancing the magnetic 
field B, 

— =Vx (vxB) + r;V2B, (6) 

ot 

In these units, v may be considered to be normalized 
by the rms value of the velocity field (C/, say), and the 
magnetic field in units that, using the square root of the 
mass density, converts a magnetic field to an Alfvcn speed 
based upon the same rms velocity. In Eqs. jSJl and lO, 



lengths are measured in units of the spherical radius R 
(the radius R will be 1 in the dimensionless units) and 
times in units oi R/U. V is the dimensionless ratio of 
pressure to mass density, with the mass density assumed 
to be spatially uniform. A forcing function f has been 
written on the right hand side of Eq. Q to represent any 
externally applied mechanical force that can be chosen to 
mimic a variety of physical effects, a common convention 
in dynamo-motivated computations in periodic boundary 
conditions. Both B and v and their curls are solenoidal. 
Note that in the incompressible case, nonlinear coupling 
between modes is provided by the nonlinearities in both 
the equation of motion and the induction equation. Here, 
we are keeping all of these. In the context of planetary 
dynamos, other physical effects sometimes justify drop- 
ping the advection term in the equation of motion. Then 
other numerical methods (see e.g. Ref. [s^ ) for linearised 
MHD can be used. 

An evolution equation for vorticity can be obtained by 
taking the curl of Eq. : 

— = V X (v X w) + V X (j X B) + i/V^w + V X f . (7) 
ot 

Given the appropriate boundary conditions, the vor- 
ticity uj can be used to determine the velocity v. The 
dimensionless viscosity v and magnetic diffusivity rj can 
be interpreted, respectively, as reciprocals of kinetic and 
magnetic Reynolds numbers based on f/, i?, and the lab- 
oratory (dimensional) values of kinematic viscosity and 
magnetic diffusivity. 

The numerical scheme is to solve Eqs. (O or {TJ and 
© by representing v and B as the Galerkin expansions, 

v(r,0 = ^^2^J,,^(r) {r<R) (8) 

qLrn 

B(r,i) = ^e,1™(i)J./m(r). (9) 

qlni 

Here, the unknowns are the now complex, time- 
dependent expansion coeflRcients and . The vor- 
ticity and current density are given by the same series, 
each term multiplied by the appropriate value of Xqi . 

Next, we may substitute Eqs. © and ^ into Eqs. 
© and 10, say, and take inner products one at a time 
with the functions J^Tm'- Using their orthonormality, 
the dynamical equations become 

j.k 

§f = E^fe^J^^-'?^'^f- (11) 

j.k 

The indices «, j, A: are regarded as shorthand; each one 
of them represents the triple of numbers q,l,m necessary 
to identify a single member of the family of the C-K func- 
tions Jqim- The sum is over all the values retained in the 
Galerkin expansion. 
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FIG. 1: Speed-up of the code in two Linux clusters in parallel 
simulations with max{g} = max{/} — 7. The dotted line 
indicates the ideal scaling. 

The nonlinearities, both in the original partial differen- 
tial equations and in the ordinary differential equations 
()10|l and are quadratic. The kinematic coupling co- 
efficients (which do not contain the ^) A^f,, Bji., and Cjj,, 
are numerical integrals of considerable complexity. Their 
evaluation and storage as a table is one of the most de- 
manding parts of the computation, and some features of 
that evaluation are described in Section llTTl On the right 
hand side of Eq. we have represented the forcing 

function f , if any, by its expansion over the retained C-K 
functions. The coefficients if non-zero, are to be re- 
garded as known, possibly time-dependent functions that 
may excite the velocity field, and may stand for various 
mechanical processes. 

The coupling coefficients in Eqs. ((TUI) and l(TT)l are 
reducible to 

^jfc — C^jk ~ ~r^k 1 Bjk — ^i^j^k y (12) 

where 

i]k^ j r^-{ijy~ik)d\. (13) 

III. NUMERICAL METHOD 

Equations l|l()() and (|ll(l are solved numerically with 
double precision in a sphere of radius R ~ 1. The ex- 
pansion coefficients ^ are in general complex, and since 
the fields are real, the coefficients satisfy the condition 
^q,i-m = (— l)™^g/m- ^ rcsult. Only the coefficients 
for non-negative values of m are stored and evolved in 
time. 

Before the simulation is started, for a given resolution 
in q and I (and all possible values of m) all required values 
of the normalization coefficients Cqi are computed using 



Eq. Q and stored. The values of \qi are computed 
numerically as the roots of the spherical Bessel functions 
using a combination of bisection and Newton-Raphson 
[33. Finally, the coupling coefficients /jj. are computed 
and stored. 

The coupling coefficients are complex, and from Eq. 
(|13|l satisfy the relation /j^, = — . The integral in Eq. 
(|13|l is separable in spherical coordinates. In the (p direc- 
tion, the integral reduces to the condition mk = rrii — rrij ; 
in any other case the coupling coefficients /j^ are zero. 
The radial integral reduces to seven integrals involving 
three spherical Bessel functions or their derivatives, and 
the integral in the polar angle reduces to seven integrals 
on three Legendre functions and their derivatives. Ra- 
dial integrals are computed numerically with high preci- 
sion using Gauss-Legendre quadratures, while inte gral s 
in 9 are computed using Gauss- Jacobi quadratures |37l| . 
Due to symmetry properties of the Legendre functions, 
all coupling coefficients with li + Ij + Ik + rui + mj + nik 
even are purely real, while the remaining coefficients are 
purely imaginary, another property that can be used to 
save memory. 

Once tables containing all these values are stored, the 
evolution of the system reduces to solving the set of or- 
dinary differential equations defined by Eqs. H1U|) and 
(|ll|l . These equations are evolved using a Runge-Kutta 
method of fourth order [s^. The MHD equations have 
three quadratic ideal invariants: the total energy E, the 
magnetic helicity H, and the cross helicity K. In spectral 
space the invariants can be computed as 

i 
i 

The conservation of these quantities up to the numerical 
precision serves as a test of the code and have been veri- 
fied in simulations with v — — 0. In a simulation with 
max{(7} = max{^} = 5, the invariants were conserved up 
to the sixth decimal place after 200 turnover times. The 
turnover time is defined as T = R/U . 

The system is evolved entirely in spectral space and all 
global quantities are also computed spectrally. To obtain 
representations of the fields in configuration space, Eqs. 
© and (jnj) are used. 

The reciprocal of the smallest |A| may be identified 
with the largest length scale in the dynamics allowed 
by the boundary conditions, and the reciprocal of the 
largest |A| retained may be considered to be the smallest 
resolvable spatial scale. In a typical computation (see e.g. 
Section|Vl), these numbers may be « 4.59 and « 41.3, re- 
spectively. The minimum and maximum wavenumbers in 
our previous 3D periodic dynamo computatior is ( e.g. in 
a 256'^ dealiased simulation using the 2/3-rule [Sg) have 
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typically been 1 and 85 in dimcnsionlcss units, by com- 
parison. The fact that the maximum to minimum ratio 
of allowed wavenumbers is more than 9 times greater for 
the 256^ code shows that far less spatial resolution ap- 
pears to be available in the C-K code as presently run. 
If "degrees of freedom" are taken as a formal measure 
of the resolution, then its 1.6 x IC* makes the C-K code 
roughly comparable with the 7.8 x 10^ available to a de- 
aliased periodic 3D code that is of resolution 64^. That 
is not, however, the whole story, since very many Fourier 
coefficients would be needed to represent accurately any 
of the C-K functions which are, individually, all consis- 
tent global physical states of the system obeying all the 
boundary conditions. We find that there are many situ- 
ations that can be computed with the C-K code that are 
nontrivial, highly nonlinear, and that are well-resolved 
as indicated by the presence of an unmistakable dissipa- 
tion range for their A-spectra; these include situations 
exhibiting a variety of dynamo behavior. 

Besides the quadratic ideal invariants and the recon- 
struction of the field components, two vector quantities 
will be of interest. The angular momentum L is defined 
as 



r X vd X , 



(17) 



where a unity mass density is assumed, and the magnetic 
dipole moment fi is given by 



(18) 



In terms of the C-K functions, these two quantities be- 
come 



fi{\Ki\R) 



3^ IVI 



+V2Re X - V2lm (Q^^O yj 

M = 2i?3y|E^«.ilVibl(|A,,il^) K^,o^- 
q 

+V2Re (e,^i,i) X - V2lm (^f y . 



(19) 



(20) 



As previously mentioned, in the code a sphere with 
i? = 1 is used. Note that only modes with / = 1 and 
m = 0, ±1 give a contribution to the angular momentum 
and the dipole moment. With the boundary conditions 
considered in this work, the angular momentum is not a 
conserved quantity unless v ~ 0. Other choices of the 
boundary conditions can lead to a conservation of the 
angular momentum even in the non-ideal case. Those 
boundary conditions apply to the case when the sphere 
of magnetofluid is isolated from torques, but we will defer 
consideration of that situation to a later paper. 

The code is written in Fortran 90 and parallelized us- 
ing MPI. Since most of the computing time is spent in the 
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FIG. 2: Total energy (dotted line), magnetic energy (solid 
line) , and kinetic energy (dashed line) as a function of time in 
run II. At late times the system is dominated by magnetic en- 
ergy. Since the only dynamic evolution possible for B requires 
a velocity and since the kinetic energy has essentially disap- 
peared, the system has "frozen" into a nearly purely magnetic 
state, which can only slowly resistively decay hereafter. 



sums in Eqs. (|10|l and (|ll|l . the parallelization is done 
as follows. Each processor contains a complete copy of 
the expansion coefficients ^, but only a portion of the 
coupling coefficients /j^. The array /j^. is distributed in 
q if the number of processors is smaller or equal than 2q, 
and distributed in q, /, m in any other case. Each pro- 
cessor computes the sums in Eqs. (fTn|) and ((TT|l locally 
for the corresponding values of q, /, m, and after each it- 
eration of the Runge-Kutta method the coefficients ^ are 
synchronized between all processors. The required com- 
munication is minimal and the scaling of the code as the 
number of processors is increased is close to ideal. 

Figure^shows the speed-up vs. the number of proces- 
sors in two different Linux clusters, in a simulation using 
max{q} = max{Z} = 7. The clusters differ in the net- 
work configuration. The speed-up is defined as the time 
required to do one time step in N processors divided by 
the time required in one processor. The code shows ideal 
scaling up to ~ 2max{g}. A drop is then observed 
and is related to the change in the parallel distribution 
of the array Jj^,. However, after this drop a linear scaling 
is again recovered. 



IV. SELECTIVE DECAY 

"Selective decay" refers to a frequently studied MHD 
turbulence situation in which a system with initial mag- 
netic helicity present evolves with a rapid decay of to- 
tal energy relative to the magnetic helicity. This situ- 
ation has been studied previously in periodic boundary 
conditions and it is of interest to see if the behavior is 
affected by spherical geometry. The late-time state is a 



6 



0.07 F 



0.06 r 



0.05 r 



0.04 



r 


/ 


; 




/ 






/I 


; 


; / 
: // 
' // 

\J/ 













1 
Time 



1 5 



20 



10" 



10" 



10"° --^ 



1 0" 




FIG. 3: Relative magnetic helicity as a function of time 
for runs I (solid line) and II (dashed line). The relative 
helicity at late times increases with the Reynolds number. 
Note the maximum possible value for the relative helicity is 
min-^{|A|} ^ 0.22. 



quasi-steady configuration in which the remaining energy 
is nearly all magnetic and is condensed into the longest 
wavelength modes allowed by the boundary conditions. 
Selective decay has been extensively studied in periodic 
boundary conditions (see, e.g., [H E Ell ) . In Ref. IH 
it was found that given isotropic initial conditions in a 
periodic box, the final state of the magnetic field corre- 
sponds to an Arn'old-Beltrami-Childress (ABC) field at 
the largest possible scale with A, B, and C equal. It is 
therefore of interest to test selective decay in a sphere, 
and to observe the geometry of the magnetic and velocity 
fields in the late-time state. 

Two runs were done, the first (run I) with v = rj = 

I X 10"^, and the second (run II) with v = rj = 6 x 10^'^. 
Run I was done with max{q} = max{l} = 7, and run 

II with max{g} = max{Z} = 9. In both simulations, no 
external force was applied, and the system was allowed 
to evolve for a long time (20 initial large scale turnover 
times). 

The non- vanishing initial expansion coefficients in both 
simulations are 



?3,3,0 — ?-3,3,0 ^ "0, ?3,3,0 ~ 3S-3,3,0 
a.3,™=C3,3,™ = "0(l + *), 
^3,3,771 = fC-3,3.m = So(l - i). 



(21) 
(22) 
(23) 



where m runs from 1 to 3 and negative values of m arc 
given by £,q,i-m = (-l)'"^,;™- The amplitudes uq and 
Bq were chosen to have initial kinetic and magnetic en- 
ergies of order unity (mq = 4 and Bq = 0.4). These 
initial conditions correspond to a small cross correla- 
tion between the velocity and magnetic fields, a non- 
helical velocity field at an intermediate scale, and a (non- 
maxinially) helical magnetic field at the same scale. The 
initial angular momentum is zero and remains negligible 



FIG. 4: Spectrum of kinetic energy as a function of A for 
run II, ai t — 0.1 (solid line), t = 1.5 (dotted line) and t = 
3 (dashed line). Note the appearance of a clear dissipation 
range. 
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FIG. 5: Spectrum of magnetic energy as a function of A for 
run II, at t = 0.1 (solid line), t — 1.5 (dotted line) and f = 3 
(dashed line). A magnetic dissipation range becomes clearly 
visible. 



during the complete simulation. The kinetic and mag- 
netic Reynolds numbers, based on the length R = 1 and 
the initial rms velocity, are respectively defined as 



Rv — 



R 



M 



RU 

V 

RU 



(24) 
(25) 



In run I Ry = Rm ~ 98, and in run II Ry = Rm ~ 165. 
Figure |21 shows the time history of the kinetic, mag- 



7 




FIG. 6: Amplitude of the coefficients ^ in run II at t = 10 as 
a function of q and / (summed over all values of m). Most of 
the magnetic energy is at the largest available scale, and an 
asymmetry is observed between positive and negative values 
of q. 



netic, and total energies in run II. At late times, the 
kinetic energy is negligible and the system is dominated 
by magnetic energy. The magnetic helicity decays slowly 
compared to the energy, as indicated by the evolution 
of the relative helicity Hm/Em (Fig- El- As time 
evolves, the ratio Hm/Em increases until a steady state 
is reached. The final state reached by the system is not a 
"Taylor state" , a state of maximum possible helicity for 
the given energy (the maximum possible value for the 
relative helicity is min~"'^{|A|} w 0.22). In run I, this 
is because after t h most of the kinetic energy has 
decayed, and as a result the spectral exchange between 
modes stopped. In run II the decay of the kinetic energy 
takes place at a later time, and as a result the final value 
of the relative helicity is larger than in run I. Larger fi- 
nal values of the relative helicity can be expected if the 
Reynolds numbers are further increased. 

Since the C-K functions are eigenfunctions of the curl 
with eigenvalue A, the Laplacian operators in the diffu- 
sion terms in Eqs. H5I7|) are proportional to A^. As a 
result, |A| plays in this case the role of the wavenumber 
k in the Fourier base. To define the energy spectrum, we 
linearly bin the spectral space in shells of constant | A| and 
sum the power of all the coefficients in each shell, in anal- 
ogy with the usual procedure in Fourier-based spectral 
methods. Figures ^ and [S] show respectively the result- 
ing kinetic and magnetic energy spectrum as a function 
of |A| in run II at three different times. 

At early times {t = 0.1), the signature of the ini- 
tial conditions in the spectrum can be easily recognized. 
Both spectra peak at |A| ~ 14, corresponding to the non- 
vanishing initial perturbation at g = 3 and 1 — 2). As time 



evolves, the amplitude of the kinetic energy spectrum de- 
cays but the position of the peak remains approximately 
constant. On the other hand, the peak in the magnetic 
energy spectrum moves to smaller values of |A|, corre- 
sponding to larger scales. At late times the system is 
dominated by magnetic energy, and most of it is concen- 
trated in the largest available scale in the domain. 

Based on the definition of the energy spectrum, we can 
also introduce an energy-containing lengthscale as 



R min{|A|} 



E 



£;(|A|)|ArM|A|, 



and a Taylor lengthscale 



At = i? min{|A|} 



Ej J i?(|A|)|Apd|A| 



1/2 



(26) 



(27) 



where E is the energy, and i^dAj) is the energy spectrum 
as a function of |A| (the sums arc represented symboli- 
cally as integrals). Using the kinetic and magnetic en- 
ergies, these characteristic lengths can be computed for 
the velocity and magnetic fields. In both runs, at i = 
£ « At ~ 0.33 for both fields. As the system evolves 
these quantities grow monotonically, but while at i = 20 
for the velocity field ^ « 0.5, for the magnetic field £ w 1. 
As a criterion to decide if the simulations were well re- 
solved in spectral space, the scales where the kinetic and 
magnetic cnstrophy spectrum peaked were observed as a 
function of time, and it was asked that their correspond- 
ing wavenumbcrs |A| were smaller than max{|A|} at all 
times. 

Figure shows the amplitude of the individual modes 
f ^ in run II at t = 10 as a function of q and I (all val- 
ues of m for each value of I are summed). Most of the 
magnetic energy is concentrated in the shell with / = 1, 
and the modes with q = ±1 in this shell have the largest 
amplitude. Note the imbalance between the mode with 
q = I and g = — 1, indicating some helicity is present in 
the magnetic field. 

The dominance of a helical large scale magnetic field 
at late times can also be identified in an inspection of 
the fields in configuration space. Figure [7| shows field 
intensity and field lines for the velocity and magnetic 
field at t = 1.5 (left) and t = 15 (right) in run II. 
While at early times both fields show small scale fea- 
tures, at late times the velocity field looks reminiscent of 
a quadrupole and the magnetic field looks like a dipole 
oriented roughly in the z direction. However, the mag- 
netic field at t = 15 is helical and the magnetic field lines 
are not purely poloidal. There is a small toroidal compo- 
nent to the magnetic field, and the magnetic field lines 
proceed slowly in the (/)-direction in a helical fashion. 

In Fig. [7| and in the following visualizations the labels 
are as follows. The x, y, and z directions are indicated 
by the arrows (in the online version, these are respec- 
tively red, green, and blue). The color and opacity are 
proportional to the field intensity, colorbars are given as 
a reference. Field lines are computed taking a snapshot 
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FIG. 7: Above: kinetic energy density and velocity field lines in run II, ai t — 1.5 (left) and at t = 15 (right). Below: magnetic 
energy density and magnetic field lines in the same run, at t — 1.5 (left) and at f = 15 (right). For convenience, intensity and 
field lines are always shown in pairs, with intensity on the left and field lines on the right. 



of the field at a fixed instant in time, and integrating a 
trajectory from twenty random points in the surround- 
ings of the center of the sphere. The field is not evolved 
in time as the lines are integrated. To indicate the direc- 
tion of the fields, in the online version the lines change 
color according to the distance integrated from the initial 
point; from red to yellow, and finally blue. 

V. DYNAMO EFFECT 

In an MHD dynamo, an initially small "seed" magnetic 
field is amplified and sustained against Ohmic dissipation 
solely by the motions of a conducting fluid. Magnetic 
fields observed in planets and stars are believed to be the 
result of a dynamo process. The mechanical mechanisms 
proposed vary widely, including thermal convection, Ek- 
man pumping due to rotation, precession, irregularities 
on the inner surface of the planetary mantles, and so on. 
In this Section, we study three simple examples of forced 
dynamo action in the sphere with relatively simple re- 
alizations of the forcing function f of Eq. (5). Figure 
|S1 shows visualizations of the two expressions used for f . 
The geometry of the forcing function is not intended to 
mimic any particular process in the planetary or stellar 
dynamos, but is rather inspired in a common practice in 
dynamo simulations with periodic boundary conditions: 
a few spectral modes are forced in the mechanical energy, 
and for Reynolds numbers large enough generic proper- 
ties of the simulations are studied. 

In all cases to be described here, a purely hydrody- 
namic (B = 0) computation was carried out first until a 
steady state (laminar or turbulent) was reached for the 




FIG. 8: Intensity of the forcing function f used in the dynamo 
simulations (left), and associated field lines (right). Above: 
function used in run III. Below: function used in runs IV and 
V. 



velocity field. The amplitude of f was chosen to have rms 
velocity of order one in the steady state. Then a random 
magnetic field with energy Em ^ 10^® was loaded into 
the modes with \q\ = I = 4. The simulations were then 
continued in order to observe the amplification and sub- 
sequent development of the magnetic field. 

Three dynamo simulations were done in which the 
Reynolds numbers and the number of modes excited were 
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FIG. 9: Kinetic (dashed line) and magnetic energy (solid line) 
as a function of time in run III. 
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FIG. 10: Trace of the magnetic dipole orientation on the unit 
sphere (above) and magnitude of the dipole moment as a func- 
tion of time (below) in run III. 

progressively increased. In the first case, the system 
reaches a steady state with a stable dipole moment which 
varies little in time. In the second case, the dipole mo- 
ment forms, then spontaneously changes direction. In 
the third, the Reynolds number is high enough that the 
velocity field might be called turbulent, and many modes 
are excited; for this case, the dipole moment changes er- 
ratically in time. 

A resolution of maxjlgj} — max{/} — 9 was used in 
all the runs. The same criteria as in the previous section 
was used to decide if a simulation was well resolved. For 




FIG. 11: Above: kinetic energy density and velocity field lines 
in run III, at t = 915. Below: magnetic energy density and 
magnetic field at the same time. The same conventions than 
in Fig. 0are used. 

the last run, a simulation with higher resolution was also 
carried out to sec if the results would be modified by the 
change in resolution. 



A. Laminar runs 

In this section we present results from two runs with 
1^ = T] = 3 X 10~^. In run III, the external forcing f in 
Eq. (O is given by the coefficient 

el,2,i =/o(l + 0, (28) 

with fo = 1.4, which corresponds to one C-K mode and 
as a result injects maximum kinetic helicity in the system. 
In run IV, the external forcing is 

^2,2,0 = 5^-2,2,0 = fo, (29) 

e2^,2,™ = 5e{2,2,,„ = /o(l+»), (30) 

where m runs from 1 to 2 and negative values of m are 
again given by ^g,;,_,„ = i~^)^^qim- ^he ampfitude of 
the forcing is fo = 0.9. This forcing injects non-maximal 
kinetic helicity (plus, of course, kinetic energy) into the 
system. Figure [S] shows visualizations of the two forc- 
ing functions in configuration space. In both cases the 
forcing is stronger in the center of the sphere, and a mod- 
ulation due to m = 2 modes in the forcing can be easily 
identified. The phase and amplitude of the external force 
f were kept constant during the entire simulations. 

Figure |51 shows the time history of the magnetic and 
kinetic energy in run III. The Reynolds numbers based 
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FIG. 12: Kinetic (dashed line) and magnetic energy (solid 
line) as a function of time in run IV. 



on the length R = 1 and the rms velocity are Rv = 
Rm ~ 290, and the energy containing scale of the flow is 
£ w 0.5. Before the magnetic field is introduced, only the 
forced mode is excited. 

After the magnetic seed is introduced, the magnetic 
energy is amplified exponentially in a kinematic regime. 
Then the magnetic field saturates around t w 150 and 
the Lorcntz force quenches the velocity field. In the final 
steady state, more mechanical modes besides the forced 
mode are excited. This is the result of an instability of 
the flow, triggered by the Lorentz force as the magnetic 
field grows exponentially. 

Figure ^| shows the trace of the orientation of the 
dipole moment in the surface of the unit sphere, and the 
magnitude of the dipole moment as a function of time 
for run III. The dipole moment grows during the kine- 
matic regime, but its orientation changes erratically in 
time. At t w 200 \fj,\ reaches its maximum amplitude 
and converges slowly to a steady state. Also it direction 
changes slowly, and after t w 400 almost no change is 
observed. At late times the dipole shows no inclination 
toward further systematic dynamical development. 

A visualization of the velocity and magnetic flelds in 
conflguration space in the steady state of run III is shown 
in Fig. 111! The kinetic energy is concentrated in two 
counter-rotating regions, located in the center of each 
hemisphere. Note the m = 2 modulation in the forcing 
is still visible in the kinetic energy density. The velocity 
field in these regions is mostly toroidal, as indicated by 
the velocity field lines. The magnetic energy is larger in 
the center of the sphere, and along the axis defined by 
the two counter-rotating eddies. In the interior, but away 
from the axis, magnetic field lines arc mostly toroidal, as 
the result of the stretching by the two counter-rotating 
eddies. Along the axis and close to the boundaries, the 
flow is mostly poloidal. 

The time history of the kinetic and magnetic energy 
in run IV is shown in Fig. E| The Reynolds num- 
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FIG. 13: Trace of the magnetic dipole orientation on the unit 
sphere (above) and magnitude of the dipole moment as a func- 
tion of time (below) in run IV. 



bcrs for this run (based on the length i? = 1) are 
Rv = Rm ~ 280, and the energy containing scale of the 
flow is I K, 0.48. Although the kinematic viscosity and 
magnetic diffusivity are the same as in run III, the ex- 
ternal forcing injects energy in a larger number of modes 
and even before the magnetic field is introduced non- 
forced modes have some mechanical energy. After the 
magnetic seed is introduced, the magnetic energy grows 
exponentially until at i « 150 saturates. The system 
seems to reach a steady state but suddenly at i sa 500 
the magnetic energy decreases by a factor of « 1.8, the 
kinetic energy increases by 1.1, and the system reaches 
a second steady state. 

The abrupt change in the kinetic and magnetic en- 
ergy in run IV at t sa 500 is associated with a reorien- 
tation of the magnetic dipole moment. Figure [T^ shows 
the time evolution of the direction and amplitude of /x. 
As in run III, in the kinematic stage the dipole moment 
grows rapidly and its direction fluctuates erratically, un- 
til reaching a first quasi-steady state aX t k, 200. The 
dipole moment then stays approximately constant until 
ai t K, 400 the magnetic field evolves rapidly, and the 
dipole moment changes direction to a second attractor 
(reached at t k, 600). The angle the dipole flips by is 
close to 7r/2. After t sa 600 only a slow change in 
is observed. The amplitude of the dipole moment also 
changes rapidly as the dipole shifts at t sa 400; while at 
t sa 300 sa 0.36, at t « 1000 « 0.54. 

Figure lTH shows the kinetic and magnetic energy spec- 
tra in run IV at different times. The kinetic energy spec- 
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FIG. 14: Kinetic energy spectrum at t = 52 [tliick (blue) solid 
line] and at t — 500 [thick (blue) dash-dotted line] in run IV; 
the thin lines correspond to the magnetic energy spectrum at 
t = 52 (solid), t = 82 (dotted), t = 250 (dashed), and t = 500 
(dash-dotted) . 

trum peaks at |A| w 9, corresponding to the forced modes 
with q = 2 and I = 2. The energy in the remaining modes 
is at least two orders of magnitude smaller than in the 
forced modes. At late times, some kinetic energy is ex- 
cited in the largest available scale, as well as a small an- 
gular momentum (|Lp/i?y « 1.5 x 10^^ after t sa 600). 
The angle between the dipolc moment and this small an- 
gular momentum remains constant after t = 100 and is 
it/ 2. Visualizations of the velocity and magnetic fields 
in configuration space arc shown in Fig. 1151 The ge- 
ometry of the velocity field is more complex than in run 
III, although the m = 2 modulation in the kinetic energy 
can still be recognized. Note a.t t = 315 the magnetic 
field lines are mostly poloidal in the center of the sphere, 
and toroidal close to the boundary. The change in the 
orientation of /x at t w 400 takes place without a strong 
change in the velocity field, and a relatively small change 
in the magnetic field configuration. 

B. A rapidly- varying dipole 

For run V, the external forcing f is given by Eqs. H29|) 
and H30|l . but the kinematic viscosity and magnetic dif- 
fusivity arc dropped to i/ = 77 = 3 x 10~*. The result- 
ing kinetic and magnetic Reynolds numbers based on the 
length R = 1 are Ry = Rm ~ 2300, and the energy- 
containing and Taylor scales are respectively I k, 0.4 and 
At = 0.35. 

The evolution of the kinetic and magnetic energy in run 
V is shown in Fig. 1161 Again, after an initial kinematic 
regime where the magnetic energy is amplified exponen- 
tially, the system reaches a statistical steady state. Note 



that in this simulation both the kinetic and magnetic 
energy fluctuate strongly with time, indicating the non- 
linear coupling between modes is stronger than in runs 
III and IV, as a result of the higher Reynolds numbers. 

Figure El shows the trace of the dipole moment in 
the surface of the unit sphere, and its intensity as a 
function of time, in run V. The direction of /x ap- 
pear to fluctuate randomly, changing hemisphere with 
a characteristic time of the order of 100 eddy turnover 
times. In the meantime, the whole surface of the unit 
sphere seems to be explored by the fluctuations in /i. 
The angular momentum is small and fluctuates around 
iLp/iJy « 5 X 10^'^. However, unlike in Run III, in this 
simulation the angle between the dipole moment and the 
angular momentum is not even approximately constant 
and fluctuates seemingly randomly between and tt. The 
question of under what circumstance the dipole favors 
one or another alignment appears to be quite unsettled, 
and deserves further study in higher-resolution simula- 
tions with larger angular momentum and rotation. 

The kinetic and magnetic energy spectra in run V at 
several times are shown in Fig. ^| More modes are 
excited in the velocity field, in accordance with the strong 
fluctuations in time observed in the kinetic energy. While 
in runs III and IV the magnetic energy spectrum peaks 
at large scales even during the kinematic regime, in run 
V at early times the magnetic energy peaks at scales 
smaller than the forcing scale. Also, after the nonlinear 
saturation of the dynamo, a fluctuation in the amplitude 
of the large scale magnetic field is observed. The minima 
arc correlated with times of minima of /x^, when the three 
components of the dipole moment fluctuate around zero. 
Most of the activity in this run is in small scales and 
fluctuations in the flow are larger than in runs III and IV. 
Even at late times when a large scale magnetic field has 
developed, intermediate scales give a large contribution 
to the magnetic energy. This also explains the strong 
fiuctuations observed in the dipole moment. Since fj. is 
proportional to the current density, the small scales give 
a large contribution to the dipole moment. 

Figure [T^ shows the magnetic and velocity field in real 
space at two different times in run V. The fields have 
more small scale structure than in run IV. Any trace of 
the TO = 2 modulation of the forcing has been completely 
lost. Also, in the kinematic regime (see e.g. the magnetic 
field at t = 75) the magnetic energy is mostly in the 
small scales, as also indicated by the magnetic energy 
spectrum. 



VI. DISCUSSION 

We have introduced some computational machinery 
that is intended for a quantitative discussion of nonlin- 
ear and incompressible magnetohydrodynamics inside a 
sphere for a wide range of initial conditions and varieties 
of mechanical forcing. It is concerned essentially with 
the analytical and computational aspects of dynamo ac- 
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FIG. 15: Above: kinetic energy density and velocity field lines in run IV, at t = 315 (left) and at t — 1065 (right). Below: 
magnetic energy density and magnetic field lines in the same run, at t = 315 (left) and a.tt = 1065 (right). The same conventions 
than in Fig. |7|are used. 




FIG. 16: Kinetic (dashed line) and magnetic energy (solid 
line) as a function of time in run V. 



tion in this geometry in a somewhat abstract setting, 
and is not the same as periodic dynamo computations 
with rectangular symmetry, or with planetary-dynamo or 
solar-dynamo computations whose central focus is repro- 
duction of observations of magnetic behavior of a real sys- 
tem. There are other features yet to be included and it is 
our intent to introduce them one at a time: rigid rotation 
and Ekman pumping, an insulating but non-conducting 
boundary to permit the magnetic field's penetration into 
the surrounding vacuum region, and different forms of 
mechanical excitations. 

The operation of the wholly spectral code, which has 
had a precedent in axially periodic circular-cylinder ge- 
ometry has been illustrated by a few simple examples. 




FIG. 17: Time evolution of the three components of the dipole 
moment in run V. Labels are as in Fig. ITol 



not in any sense a comprehensive study. First, decay- 
ing turbulence has been studied involving the relaxation 
of helical initial conditions to a magnetically dominated 
state whose spectrum is dominated by the longest al- 
lowed spatial scales and whose kinetic energy is essen- 
tially gone. Also, dynamo computations have been done 
for three successively higher sets of Reynolds numbers, 
revealing a different magnetic behavior in each case. In 



13 



1 .000 




10 

X 



FIG. 18: Kinetic energy spectrum at t = 75 [tliick (blue) solid 
line] and ait — 600 [thick (blue) dash-triple dotted line] in run 
V; the thin lines correspond to the magnetic energy spectrum 
at i = 75 (solid), t = 168 (dotted), t = 240 (dashed), t = 544 
(dash-dotted), and t = 500 (dash-triple dotted). 

the first case, a magnetic dipole formed in what was es- 
sentially a laminar velocity field and appears willing to 
persist for as long as the code is run. In the second 
case, another dipole formed, but after the passage of a 
few hundred large-scale eddy-turnover times, it changed 
its magnitude to a somewhat larger value, and flipped 
its orientation, for reasons we do not know but that are 
worth exploring. These fluctuations have required nei- 
ther a preferred direction enforced by rigid rotation nor 
thermally convective rolls. Finally, in the third case, the 
velocity field had Reynolds numbers (based on the ra- 
dius of the sphere) of about 2300 and could be said to 
be turbulent. In this turbulent case, there was a dipole 
moment, but it exhibited no systematic or regular be- 
havior as far as we could tell, and changed its orientation 
from one hemisphere to the other every 100 or so eddy 
turnover times, exploring in the meantime all possible 



orientations. It did appear to be a resolved computa- 
tion, despite the turbulence, and we believe represents 
a bona fide solution of the MHD equations. All these 
runs, it should be stressed, were carried out with a mag- 
netic Prandtl number of unity, and may well change for 
values far from that. One outcome to be noted is that 
neither turbulence nor rotation have been a necessary in- 
gredient for the development and maintenance of a mag- 
netic dipole, but the presence of mechanically helicity has 
helped a lot. Spontaneous changes in magnetic dipole ori- 
entation have been easy to observe, both in laminar and 
turbulent cases. 

We should also stress that the code as presently consti- 
tuted is limited to a relatively small number of degrees of 
freedom, far fewer than wcU-rcsolvcd high Reynolds num- 
ber mechanical turbulence would demand. It will also be 
attempted to design wall-friction terms to permit a more 
efficient exchange of angular momentum from the fluid to 
the wall |30j| . once rigid rotation is introduced. Also rigid 
rotation and different boundary conditions for the mag- 
netic fields can be implemented. Given the properties 
and limitations of the method described (purely spec- 
tral, non-dispersive, and conservative), we believe this 
method can be used as a testbed to explore the effect of 
different physical effects, boundary conditions, subgrid 
models (several models, such as the Lagrangian averaged 
model 0, 0| are easier to implement in spectral 
space), etc., before trying these ideas in more complex 
and realistic codes to reach high Reynolds numbers. 
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